function [r, p, z, df] = pearson(x, y)
% compute peason's product moment correlation, the p value, and
% a usefule z index from x and y for a linear correlation between x and y
% from p638 Press et al.
%
r=[];
p=[];
z=[];
n=length(x);
m=length(y);
if(n ~= m)
	disp('Lengths of input vectors must match')
	return;
end

if(n < 1)
	disp('Cannot fit line less than 2 points')
	return;
end;

TINY=10^-20;

ax=mean(x); % compute means of the variables
ay=mean(y);

xt=x-ax; % difference from mean value
yt=y-ay;

sxx=sum(xt.^2); % sums of squares etc.
syy=sum(yt.^2);
sxy=sum(xt.*yt);

r = sxy/sqrt(sxx*syy);
z=0.5*log((1.0+r+TINY)/(1.0-r+TINY)); % fisher's z transform
df=n-2;
t=r*sqrt(df/((1-r+TINY)*(1+r+TINY)));
%p=beta_i(0.5*df, 0.5, df/(df+(t*t)));
p = erfcc(abs((z)*sqrt(n-3))/sqrt(2));
return;
